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Abstract 

We introduce a ra-term quadrature to integrate inner products of n functions, as 
opposed to a Gaussian quadrature to integrate 2n functions. We will characterize 
and provide computataional tools to construct the inner product quadrature, and 
establish its connection to the Gaussian quadrature. 
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1 The inner product quadrature 

We consider three types of n-term Gaussian quadratures in this paper, 

• Type-1: to integrate 2n functions in interval [a, b\. 

• Type-2: to integrate n 2 inner products of n functions. 

• Type-3: to integrate n functions against n weights. 

For these quadratures, the weight functions u are not required positive definite. Type-1 
is the classical Guassian quadrature, Type-2 is the inner product quadrature, and Type-3 
finds applications in imaging and sensing, and discretization of integral equations. 

In this section we will introduce and characterize the Type-2, inner product quadra- 
ture. $2]presents algorithms for constructing the quadrature, ^establishes framework to 
link the first two types, and introduce the Type-3 quadrature. §H illustrates our quadra- 
ture design method with several examples. $5] explores generalizations of our quadrature 
methods to higher dimensions and examines applications to inverse scattering problems. 



1.1 Notation 

For x G [a, b] and k G [a, by the usual abuse of notation we will denote by G(k,x) 
the four related objects 

1. A family of L 2 functions on [a, b], with k G [a, (3] 

2. The linear subspace spanned by these functions; 

3. The kernel of an integral operator; 

4. The matrix of that operator of size [a, f3] -by- [a, b\. 

When k in G(k,x) takes on some finite n values kj in the resulting set of n 

functions are denoted by 

T(n,x) = { G(k,x), x G [0, 1], k = kj, j = 1 : n} (1) 
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Viewed as a matrix of size n-by-[a, b], T(n, x) has infinite number of columns, and its n 
rows consist of the n functions Tj(x) in L 2 [a,b\. For example, when k ranges from to 
n — 1, the power functions 

G(k,x) = {x k , are [0,1], k E [a, /?] } (2) 

becomes T(n, x) = IT n , polynomials of degree less than n. Similarly, when k takes on 
integers, the pure tones 

G(k,x) = { exp(ikx), ze[-7r,7r], fee [-/?,/?]} (3) 

becomes trigonometric polynomials. We will first consider a n-term quadrature to inte- 
grate the inner products in T(n,x). 

For c e [a, b], the column of T(n, x) taken at x = c is the n-by-1 vector 

T(n, c) = T(n,x)\ x=c (4) 

Given a weight function u we define the dot product in T(n, x) 

f ■ g = u(x)f(x)g(x)dx (5) 



When u is positive, (jSJ) will be adopted as the inner product for L 2 [a, b\. 

Given n distinct points Xj, j — 1 : n in [a, b], let T(n, {xj}) be the n-by-n matrix 
formed by the n columns T(n,Xj), j = l:n. 

Let B(n,n) be the n-by-n Gramian matrix of the n functions Tj(x) so that 

B(n,n) = T(n,x) ■ T(x,ri) (6) 

where T(x,n) is the complex transpose oiT(n,x). Note that the dot sign requires inner 
product (JS} with the underlying weight function u. For real valued u, B is Hermitian. 
For positive u, let T + (x,n) be the pseudo inverse of T(n,x) : L 2 [a,b] i— > C n so that 

T + (x,n) = T(x,n)B~ 1 (n,n); (7) 

thus P n : L 2 [a, 6] H- L 2 [a, 6] defined by the formula 

P n (x,y) = T + (x,n)T(n,|/) (8) 

is the orthogonal projector onto T(n, x) 
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1.2 A n-term quadrature for inner products 

For a positive u, let Q(n,x) denote n orthogonal basis functions for T(n,x). 

Definition 1.1 A n-term quadrature {xj, Wj} is one with distinct Xj in [a, b] and nonzero 
Wj, j = l:n. 

Theorem 1.2 ("Duality of row and column orthogonalities^) Let the weight u 
of ([3J) be positive. There is a n-term quadrature {xj,Wj} to integrate all inner products 
in T(n, x) if and only if the n columns of the n-by-n matrix Q(n, {xj}) are orthogonal. 

Proof Obviously, the n-term quadrature {xj,Wj}, if exists, integrates the Gramian 
Q(n,x) ■ Q(x,n) = I with positive weights Wj, namely 

I = Q(n, {x i })diag{w j }Q({x j }, n) (9) 

so the n-by-n matrix Q(n, {x J })diag{ v /uJj} is unitary and thus the columns of Q(n, {xj}) 
are orthogonal. 

Now assume that the columns of Q(n, {xj}) are orthogonal. Let the norm of the j-th 
column be 1/ y/Wj so that Q(n, {x J })diag{y / M7J} is unitary, which implies that (Q holds, 
namely there is a n-term quadrature {xj, Wj} to integrate the Gramian Q(n, x)-Q(x, n) = 
I; therefore it integrates all inner products in T(n,x). 

2 Construct the inner product quadrature 

It is a difficult, nonlinear problem to select n orthogonal columns out of infinite number 
of columns of matrix Q. The selection process can be made a lot easier if the n-term 
quadrature is requried to do a bit more. In addition to B, if the quadrature also integrates 
another Gramian 

A(n,ri) = T(n,x) ■ n(x)T(x,ri) (10) 

where \i is a simple function, then the quadrature nodes Xj will be recorded in \i. In 
fact, fi(xj) will be eigenvalues of the quotient matrix AB~ l . We first formulate this fact 
in §2. II for polynomial T(x,n). The general case is treated in §2.21 

2.1 Polynomial case 

Let T(n, x) be the n dimensional space Il n for polynomials of degree less than n. Let 
fi(x) = x so that A(n, n) is given by 

A(n,n) = T(n,x) ■ xT(x,n) (11) 
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Theorem 2.1 If there is an-term quadrature {xj,Wj} to integrate the Gramians A and 
B, then the nodes Xj are eigenvalues of AB^ 1 

Xj = \j(AB~ l ), j = l:n (12) 

provided that B is invertible (when u is not positive definite). 

Proof The n-term quadrature exact for B is of the form 

B = T{n, {x J })diag{w i }T({x i }, n) (13) 

It follows from invertibility of B that the quadrature weights w have no vanishing entry 
and T(n, {xj}) is invertible; therefore, 

AB- 1 = [T(n, {x j })diag{ Wj x i }T({x j },n)][T(n, {x j })diag{ M ; j }T({x j }, n)]" 1 (14) 
= Tin^x^dmgix^Tin^x,})}- 1 

Thus, the j-th eigenvalue of AB" 1 is Xj with the eigenvector T(n, {xj}). 

If the weight function u is positive, then by the proof of Theorem 11.21 the quadrature 
weights are 

Wj = \\Q(n, Xj )\\l (15) 
If u is not positive definite, since the quadrature is exact for the first column of B, Wj 
will be determined by solution of the n linear equations for w 

T(n, {x i })diag{t/7 i }T({x i }, 1) = T(n, x) ■ T{x, 1) (16) 

where T({xj}, 1) is the first column of T({xj}, n), and T(x, 1) = T x (x) is the first column 
of T(x, n). 

Theorem 2.2 Let the weight function u be positive definite, and let Vj denote the j-th 
eigenvector of a matrix. There is a n-term quadrature {xj, Wj} to integrate the Gramian 
matrices A and B if and only if 

XAAB' 1 ) = xj, j = l:n (17) 
Vj{AB~ l ) = T(n, Xj ), j = l:n (18) 

Proof Only need to consider orthonormal basis T(n, x) for which AB -1 = A is Her- 
mitian with orthogonal eigenvectors. Hence the proofs of Theorems 12.11 and 11.21 can be 
adopted to establish necessity and sufficiency of (jT7|) . ( fT8|) 

Obviously, the n-term quadrature {xj,Wj} integrating the two Gramians integrates all 
polynomials of degree less than 2n. Thus, Type-1 and Type-2 quadratures are the same 
for the polynomial case. 
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2.2 Arbitrary functions 



In this section we will characterize and construct an inner product quadrature for a set 
of arbitrary functions T(n, x). 

Theorem 2.3 Let B be invertible, and let Xj,Vj denote the j-th eigenvalue and vector of 
a matrix. If there is a n-term quadrature {xj,Wj} to integrate the Gramians A of ( TJ7| ) 
and B, then 

X^AB- 1 ) = n(xj), j = l:n (19) 
v j (AB~ 1 ) = T(n,Xj), j = l:n (20) 

The proof is identical to that of Theorem 12.11 

Definition 2.4 A function \i is said to be a minimal function ofT(n,x) if 

r(T,/x) =: rank{[(J - P n )fx(x)T{x,n)}} = 1 (21) 

In other words, the n functions fi(x)T(x,n) don't entirely lie in the span of T(x,n), but 
the part of n(x)T(x, n) that is outside of T(x, n) is required to be minimal - the residual 
dimension r(T, /i) is 1. For example, \x = exp(ix) is a minimal function for the subspace 

E m = span[exp(i/cx), k = —m : m], m > 0, x G [— 7r,7r] (22) 

More general definition for minimal function will be given in §3.21 Modifications are also 
necessary for higher dimensions. 

Theorem 2.5 Let the weight function u be positive definite, and let Vj denote the j-th 
eigenvector of a matrix. The three conditions are equivalent 

(i) There is a n-term quadrature {xj,Wj} to integrate A and B 

(ii) The quotient matrix AB -1 is diagonalizable with 

\ 3 {AB~ l ) = M(^), j = l:n (23) 
v 3 (AB~ l ) = T(n,x,), j = l:n (24) 

(Hi) There exist such {xj} that for every p n G span[(J — P n )n(x)T(n, x)} 

p n {xj) = 0, j = l:n (25) 

and that the n-by-n matrix T(n, {xj}) is invertible. 
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Proof The proof of equivalency of (i) and (ii) is similar to that of Theorem 12.21 Now 
we establish equivalency of (ii) and (iii). By (ii), and by (J7|) and (jSJ), 

n(xj)T{n,Xj) = A_B _1 T(n, Xj), j = l:n (26) 
= T(n, x) ■ fi(x)T + (x, n)T(n, xj)) 
= [T(n, x)fi(x)] ■ P n (x, xj) 
= {P n [T(n,x)fi(x)]} x=Xj 
= {(I- I + P n )[T{n,x)n{x)]} x=Xj 
= T(n,Xj)n(xj) - {(I - P n )[T(n,x)fi(x)]} x=Xj 

which holds if and only if 

{(I-P n )\pL(x)T(n,x)]}( Xj )=0, j = l:n, (27) 
namely ( 1251) holds, and the n-by-n matrix T(n, {xj}) is invertible. 

Theorem 12.51 does not require that fi is minimal, but if it is not then all p n G span[(J — 
P n )fj,(x)T(n, x)] must share n common roots at the quadrature nodes Xj. In other words, 
it is unlikely for a n-term quadrature to integrate both B and A exactly if /i is not 
minimal. 

If T(n,x) = Il n , polynomials of degree less than n, then /i(x) = ax + f3, a ^ is 
a minimal function whereas x 2 is not. With a minimal fi, p n of (12 5 j) is the orthogonal 
polynomial of degree n, provided that the weight u is positive definite. The condition 
( 125]) is well known as a part of the Gauss formula. 

However, when the functions T(n, x) are not polynomials, the n-term quadrature for- 
mula may not be a Gaussian quadrature in the classical sense. In general, to integrate the 
inner products in A and B is not the same as to integrate some 2n functions. Conversely, 
given a set of 2n functions to integrate by a Gaussian quadrature, additional work is re- 
quired to reformulate this Type-1 quadrature as a Type-2, inner product quadrature. 
This issue will be addressed in the next section. 

3 Product law and minimal functions 

In this section we will establish framework for converting the Type-1 quadrature to 
Type-2, inner product quadrature. While an inner product quadrature is natural in 
its own right and immediately useful in many applications, other applications require 
Type-1 quadratures. For many familiar and widely used families of functions the two 
quadrature problems turn out to be equivalent or nearly so. We introduce the notion of 
factor space in §3.1l and minimal function in §3.2l to connect the two types of quadratures. 
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3.1 Factor spaces 

In this section we introduce the product law and factor space for a given set of 2n 
functions, so as to convert a Type-1 quadrature for the 2n functions to a Type-2 for the 
Gramian matrix of the factor space. 

Let the rows of V(2n, x) : L 2 [a,b] t-» C 2n consist of a set of 2n linearly independent 
functions, which span a linear subspace of L 2 [a, b] denoted also by V(2n,x). 

Definition 3.1 The linear subspace V(2n, x) is said to have a factor space T(n, x) with 
a multiplier fi if 

span{Ti(x)T,,(x), T^/i^T^x), 1 < i, j <n} = V(2n,x) (28) 

As an example, the linear space H.2n for polynomials of degree less than 2n has a factor 
space Il n with fi(x) the multiplier. Likewise, let 

G m = spanfl, sm(jx), cos(jx), j = l:m — 1], m > 1 (29) 

Then G m is a factor space of G^m with /i(x) = cos(x), whereas E m of ( 122|) is a factor 
space of E 2 m with /i = 1. 

Obviously, a quadrature integrating the inner products in the factor space will also 
integrate the functions in V(2n,x). In this respect, the notion of a factor space can be 
relaxed in two directions (i) Let the span in (12 8j) include, rather than equal to, V(2n,x) 
(ii) Let the span in ( 1281) approximate V(2n, x) to a given precision. 

Definition 3.2 The linear subspace V(2n, x) is said to obey the product law if there exist 
n functions T(n,x) such that V(2n,x) is a subspace of the product space 

n(T) = spanjT, (£)!}(£), 1 < i,j < n) (30) 

A linear subspace V(2n,x) is said to obey the product law to precision e > if for any 
f G V(2n, x) the distance between f and Tl(T) is e. 

As an example, by Neumann's addition formula 9.1.78 of [T], 

m oo 

J m {x) = Jk{x/2)J m _ k {x/2) + 2 ^(-l) fe J k (x/2) J m+k (x/2) (31) 

fc=0 fc=l 

For x G [0, b], and for a prescribed precision e > 0, there exists 5 > so that 

|J s (x/2)| < e, s>b + 5 (32) 
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Thus, to precision 0(e), only finite number of terms in ( )3~Tj) remain: J s (x/2) for < s < 
b + 5. Consequently, the space 

V = span{ J m {x), < m < 2(6 + 5)} (33) 

obeys the product law to precision 0(e), and a quadrature integrating the inner products 
in 

T = span{ J,(x/2), < s < b + 6} (34) 

exactly or to precision 0(e) will also integrate functions in V to precision 0(e). 

Factor space and product law can also be extended to a family of infinite number of 
functions, denoted by G(k,x), x G [a, b], with k G [a, /3] the family parameter. 

Definition 3.3 T/ie family of functions G(k, x) obeys the product law if there exist two 
families of functions T(k,x), S(k,x), x G [a, b], k G [ai,/3i], k G [a^,/^] such that 
G(k, x) is a subset of the product space 

U(S,T) = sp&n{T(k,x)S(n,x), k G [ai,/3i], k G [a^,/^]} (35) 

Moreover, G(k, x) is said to have factor spaces T(k, x) and S(k, x) if 

span{G(A;, x), a < k < 0} = U(S, T) (36) 

Finally, G(k, x) is said to have a factor space T(k, x) if it has the factor spaces T(k, x) 
and S(k, x) with S(k, x) = T(k, x) . 

Accordingly, the n-by-n Gramians B of and A of ffTUl) can be extended to operator 
case. 

Definition 3.4 Given a function n(x), the linear operators defined by 

A(k,k') = T{k,x) ■ /i(x)S(x,k') (37) 
B(k, k') = T(k,x)-S(x,k') (38) 

are referred to as the Gramians associated with the factor spaces T(k,x) and S(k,x). 

Theorem 3.5 Let the m-by-n matrices 

A(m,n) = T(m,x) ■ fi(x)S(x,n) (39) 
B(m,n) = T(m,x) ■ S(x,n) (40) 

be the Gramians associated with the m-by-[a, b] matrix T(m, x) and the n-by-[a, b] matrix 
S(n,x) and a scalar function ji. Let the rank of B be r. If there is a r-term quadrature 
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{xj, Wj j = l:r} precise for A and B, then the m-by-m matrix AB + has r eigenvalues 
and corresponding eigenvectors of the form 

Xj(AB+) = ufa), j = l:r (41) 
vj(AB + ) = T(m,Xj), j = l:r (42) 

The remaining m — r eigenvalues are zero. 

Proof The proof is similar to that of Theorem 12.11 Since B is of rank r, the existence 
of the r-term quadrature, exact for B, implies that the quadrature weights w have 
no vanishing entry and the m-by-r matrix T r = T(m, {xj}) and the r-by-n matrix 
S r = S({xj},n) are both full rank; therefore, 

B + = [T r diag{ Wi } S r ] + = S+ diag{ Wj j' 1 T+ (43) 

Using S r Sr = I we have 

AB + = T r {n(x j )}diag{w j }S r S+diag{w j }- 1 T+ (44) 
= T r {/i(x j )}T+ 

It follows immediately that the m-by-m square matrix AB + , being of rank r or less, will 
have m — r zero eigenvalues, and owing to T r + T r = / the r remaining eigenvalues and 
vectors are given by ( 14T]) . (]42p . 

Theorem 12. II is a special case of Theorem 12 . 3 1 which is a special case of Theorem 13.51 For 
quadrature design, we are only interested in eigenvectors, if exist, of the form T(m, Xj). 

Definition 3.6 The eigenvectors of AB + of the form T(m,Xj) are referred to as the 
position eigenvectors. 

The existence of the position eigenvectors is necessary for that of a Gaussian quadrature. 
The next theorem, straightforward to verify, says that the eigenvalues for the quotient 
matrix is invariant under the change of bases by ( 1451) , ( 1461) . 

Theorem 3.7 Let the square matrices t(m,m) and s(n,n) be invertible. Let the change 
of bases, from T(m,x) to T(m,x), and from S(n,x) to S(n,x), be defined by 

T(m,x) = t(m,m)T(m,x) (45) 
S(n,x) = s(n,n)S(n,x) (46) 

Then the two quotient matrices AB + associated with the old and new bases are similar, 
with t(m, m) as the similarity transform. 
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3.2 Minimal functions 

Minimal function was defined in §2.21 for a set of n functions T(n, x). In this section, we 
will introduce minimal function for a family of infinite number of functions G(k, x). 

Definition 3.8 (Informal) A method to grow a family of functions G is to multiply the 
existing family members by a function fi, which may not be in the family. The resulting 
functions are linearly combined with those in the family to generate a new function. The 
function \i, with proper normalization, is the minimal function. 

Typical 3-term recursions use this scheme to generate a class of functions. For example, 
the Bessel functions require [i[x) = 1/x as the multiplier to push the family one step 
forward, or backward. 

For a precise definition of minimal function, let the new function G(f3 + h, x) be 
generated by linear combination of fi(x)G(k,x) and G(k,x) over k G [a, j3]. We scale \i 
such that it appears in the linear combination as follows 

G(/3 + h,x) = hG(/3, x)n(x) + G(/3, x) + tail (47) 

The tail vanishes as h — > 0, provided that /i is the log derivative of G with respect to k. 

Definition 3.9 Let G(k,x) be differentiate with respect to k in [a, 0] for almost every 
x G [a, b] . The function 

fj.(x,k)\ k=p = |^logG(fc,a:)| (48) 

is referred to as a specific minimal function ofG(k,x) at k = (3. If /i(x, k) is independent 
of k or if the dependence is separable 

/i(x,k) — p(k)q(x) so that dkG(k,x) = p{k) G(k,x) q(x) (49) 

then it is referred to as the (general) minimal function ofG(k,x). 

By (I48p . the minimal functions for the power functions ([2]) and exponentials ([3]) are 
log(x) and x. By /i's dependence on k, we divide G into three varieties 

(V.l) It is independent of k. 

(V.2) The dependence is separable. 

(V.3) The dependence is not separable. 

There are two cases for constructing a quadrature, whether Type-1 or 2 



11 



(C.l) Design a quadrature with a given weight u. 



(C.2) Design a quadrature without u given explicitly. 

(C.l) is typical of quadrature design for numerical integration; the weight u is given 
explicitly. (C.2) arises from certain applications such as inverse problems or signal pro- 
cessing where the measurement or signal is the exact integrals 

s(k) = j G(k,x)u(x)dx, ke[a,(3] (50) 

J a 

with an underlying, fixed, but unknown u. 

For (C.2), the only data available for Type-1 quadrature design is s(k). When re- 
formulated as a Type-2 quadrature problem, (V.l) and (V.2), not (V.3), will be useful 
in constructing the Gramians A and B out of the data s(k). The procedures for con- 
structing the Gramians by (V.l) and (V.2) are so similar that in the sequel we will only 
consider (V.l), namely (V.2) with p{k) = 1. 

For (C.l), the weight function u is given and the Gramians can be constructed directly 
by their definitions (137|) and (138]) for a Type-2 quadrature, or for a Type-1 quadrature 
provided that G has factor spaces T(k,x) and S(k,x). (V.3) will be useful for (C.l). 

A specific minimal function exists for an arbitrary system of functions G(k,x). In 
contrast, only certain function classes have (general) minimal functions. The next theo- 
rem is a direct consequence of Definition 13.91 

Theorem 3.10 G(k,x) has a minimal function if and only if 

G(k,x) = exp(p(k)q(x))r(x) (51) 

Moreover, ifG(k,x) has a minimal function then it has a factor space 



T{k,x) = [G(k,x)} 1/2 = ex V {p{k)q{x)/2)^r~[x) (52) 

For example, the family x k = exp(k log(x)) is of this exponential type. The family 
k x = exp(log(/c)x) is also of this type. 

When k takes only on discrete values, say integers, the differential form ( |48|) for fi 
can be replaced by a finite difference for certain classes of functions, among them are 
polynomials and modified Bessel functions: 

u(x) = = X — 1 (53) 

ry ' 1 • x n 17 
M*) = ^-7^ = ~n/x (54) 
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where x — 1 can be normalized to x, and the separable /i(x,n) = —n/x to 1/x. 

For certain applications k G [a, /3] is restricted on a uniform grid with step size h, so 
that only a finite number of the family members G(a + jh, x), j — : n, nh — (/3 — a), 
are to be integrated. For example, 

= x h , for power functions (55) 
nix) = e lhx , for exponentials (56) 
n(x) = cos(x), for trigonometries (57) 

are appropriate minimal functions. 

Remark 3.11 The minimal function introduced in Definition \3.9\ is asscociated with 
differentiation with repect to k. In general, in ( TJP| ) is replaced by a map L k which op- 
erates on G(k,x) as a family of functions ofk. For example, if Sturm- Liouville equation 
{—L x + k 2 )u(x) = has a solution of the form u = G(k, x) = G(kx), then interchanging 
the roles of k and x we have LkG(kx) = x 2 G(kx). In other words, the family of functions 
G(k,x) = G(kx) has a minimal function x 2 with respect to the operator Lk- 



3.3 Fold data into Gramians - signal processing 

In quadrature design for numerical integration, the weight function u is usually pre- 
scribed. For other applications, such as optimal design or inverse problems, u is either a 
variable or not given explicitly. 

When u is not given and the exact integrals s{k) of ( 150]) is the only available data, 
the first step in quadrature design for G(k,x) is to process the signal s to construct the 
Gramians A and B. 

In this section, we will describe the signal processing operations for converting Type-1 
quadrature for G to Type-2 for the Gramians. This signal processing is not required to 
contruct the Gramians if u is available. 

Let G(k, x) have a factor space T(k, x), k e [ai, 0i] and a minimal function \i so that 

span{G'(fc, x), a < k < (3} = span{T(A;, x)f(k', x), a x < k, k' < ft} (58) 

and 

s'(k) = / n(x)G(k,x)u(x)dx (59) 



with the latter obtained by differentiating ( 150]) . By ( 158]) . there exists linear combination 
coefficients F(k, k', n) to reproduce T(k, x)T(k', x) as a linear combination of G(k, x)\ 

T{k,x)T{k',x) = F{k,k',K)G{K,x)dK (60) 

J a 
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Integrating (1601 with respect to x against fi(x)u(x) we rewrite the result and (|59|) in 
matrix form 

s'(Jfe) = G(fc, x)fi(x)u(x) (61) 

T(k,x) ■ fi(x)T(x,k') = F(k,k',K)G(K,x)n(x)u(x) = F(k,k', k)s'(k) (62) 

The operator T(fc, x) • fi(x)T(x, k'), by ( l37l) . is Gramian matrix A{k, k'), and ( |62l) shows 
that the derivative of the signal is required to construct A, and that how the vector 
s' is packed into A by the folding operator F. The other Gramian matrix B is also 
constructed by the same folding process performed on the signal s 

A{k,k') = F(k,k',K)s'(K) (63) 
B(k, k') = F(k,k',K)s(K) (64) 

As an example, let G{k,x) be exponentials defined by ([3]), which has a factor space 

T(k, x) = { exp(ikx), x E [-it, tt], ke [-0/2, (3/2} } (65) 

By ( |60l) . the folding operator is 

F(k,k',K) = 5(k-k' -«), fc,fc' G [-P/2,/3/2], k&[-P,P\} (66) 

so that for a fixed k, the kernel -F(/c, fc', k) is zero everywhere except on the diagonal 
k — k' = k; the Gramians A and 5 of ( 1631) and (1641) are Toeplitz matrices with s'(«) and 
s(k) on the diagonal k — k' = k. 

Folding a data vector, or signal, into a matrix or matrices and subsequently processing 
them is inherently a data analysis procedure. When G(k,x) and its factor space T(k,x) 
share the same minimal function /x, the Gramian matrices A and B are of the form 

B(k, k') = T(k,x) -T(x,k') (67) 
A(k, k') = T(k,x) ■ fi(x)T(x,k') (68) 
= [d k T(k,x)}-T(x,k') = d k B(k,k') (69) 



Theorem 3.12 Suppose that G(k,x) and its factor space T(k,x) share the same mini- 
mal function fi, and that the weight function u of |3P is nonzero almost everywhere on 
[a, b}. Then the quotient matrix Q = AB~ l is the differential operator, with respect to k, 
restricted on the subspace T{k,x). 

If u vanishes on a subset of [a, b] of positive measure, the quotient matrix Q will still be 
a differential operator restricted on the range space of B, which is a subspace of T(k, x). 
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3.4 Regular izat ion 



Once the two Gramians are constructed, there are two issues with computing the quotient 
matrix Q = AB^ 1 (i) Inverting the compact operator B (ii) For a prescribed precision 
e > 0, replace Q by a finite, n-by-n square matrix for subsequent eigen decomposition. 
The two issues can be tackled together by regularization of A, B: Approximate A, B 
with finite rank operators A n , B n . 

Ideally, we should find a function s n (k) to approximate the data s(k) in a least squares 
sense to the prescribed precision which when packed by ( 1641) gives rise to B n of rank n. 
Solving such a nonlinear least squares problem is not known to be tractable in cost or 
convergence, so suboptimal schemes are sought instead. One of them requires SVD on 
B with e as the cut off precision to construct a rank n best approximation to B, so that 

B(k, k') « U(k, ra)E(ra, n)V(n, k') (70) 
A(k, k') « U(k, n)S(n, n)V(n, k') (71) 
Q(k,k') « U (k, n)S(n, n)E _1 (n, n)U(n, k') (72) 

namely, both B and A are projected on the n dimensional column (or range) subspace 
spanned by U(k,n) and row (or domain) space spanned by V(k,n). Note that while E 
is diagonal, S is generally not. Finally, by Theorem 12.31 a ra-term quadrature of finite 
precision proportional to the prescribed can be attempted by solving the eigenvalue 
problem for the projected version of Q 

Q(n, n) = S(n, n)£ -1 (n, n) (73) 



3.5 Type-3 quadratures for integral equations 

Let A, B be the m-by-n Gramian matrices of Theorem 13.51 Let the rank of B be r. A 
r-term Type-3 quadrature uses the nodes {xj,j — 1 : r} and weights W = {wij % — 1 : 
m,j = 1 :r} to integrate A and B 

A(m,ri) = W(m, j)fi(xj)S(xj,n) (74) 
B(m,n) = W(m, j)S(xj,n) (75) 

In other words, the m functions T(m, x) of (1371) are regarded as the weight functions for 
the Type-3 quadrature. 

Theorem 3.13 If there is a r-term quadrature f/4\ ), (T75| ), then them-by-m matrix AB + 
has r eigenvalues and corresponding eigenvectors of the form 

\j{AB + ) = n{ Xj ), j = l:r (76) 
Vj(AB + ) oc W(m,j), j = l:r (77) 
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The remaining m — r eigenvalues are zero. 

The proof is nearly identical to that of Theorem I3.5[ and is omitted. Let 

v(y)= G(y,x)u(x)dx, ye[a,b] (78) 

J a 

be an integral equation for u on [a, b]. Let {yi,i = 1 : m} be m points in [a, b}. Let 
T(m, x) = {G(jji, x), i = 1 :m}. Finally, let u G S^x, n), namely u is in the span of the n 
functions 5*. Then W of Theorem 13.131 is a discretization of the integral equation 

r 

v( yi )^Y, W ^) (79) 

which is precise for u G S(x,n). 



4 Examples 

In this section we present several examples to illustrate our quadrature design meth- 
ods. In §4.11 we construct quadratures for non-positive definite weight u. §4.21 and §4.31 
construct quadratures for power and exponential functions. 



4.1 Quadratures for non-positive definite weights 

Gaussian quadratures may not exist for non-positive definite weights. As an example, 
we consider n-term Gaussian quadratures to integrate polynomials of degree less than 
2n, against the weight function 

u(x) = sin(37rx) (80) 

in [—1, 1]. The oddness of u and the optimality of Gaussian quadrature preclude Gaussian 
quadratures of odd n, otherwise x = must be a quadrature node where u vanishes which 
makes the node useless. Not all even n values support a Gaussian quadrature. For the 
weight given by (I80p . there is a Gaussian quadrature for n = 16, and n = 18, but not for 
n = 14. Whenever there is a Gaussian quadrature, it can be constructed by Theorem 
12.11 Figure [1] shows the locations of the quadrature nodes in [—1, 1], and the quadrature 
weights. The weights are negative wherever u is negative. 

4.2 Power functions, Hankel Gramians 

To integrate the power functions 

G(k,x) = {x k , xe[0,l), ke[a,0\} (81) 
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Quadrature nodes for n=18, weight function sin(37i x) 




Figure 1: Quadrature nodes and weights for u(x) = sin(37rx) 



against a weight u with a n term quadrature, we follow §3.31 to construct Gramian B 
from the exact integrals s(k), and Gramian A from s'(k). 

By (|48p . the minimal function is /a(x) = log(x). The power functions obey the product 
law of Definition 13.31 with 

T{k,x) = {x\ xG[0,l], k G [a/2, 0/2] } (82) 

By flU, the folding kernel F, cf (j66]h is 

F(k,k',K) = 5(k + k' -k), k,k' G [a/2,/3/2], «e[a,/?]} (83) 

Therefore, the Gramians A and i? are Hankel matrices with s'(k) and s(k) on their 
anti-diagonals fc + k! = k. 

For a numerical experiment, we construct a Gaussian quadrature for G(k,x) = x k , 
x G [a, b] = [0, 1], /c G [a, /3] = [—1/3, 1/2] by constructing an inner product quadrature 
for the factor space T(k,x) = k x , x G [a, b] = [-3,3], k G [a/2, ,5/2] = [-1/6,1/4]. 
Following the procedures of §3.41 a n term quadrature, though not precise to integrate 
all functions in G(k, x), was constructed from the n-by-n Gramians A and B of (1701 For 
a cut off precision e = 10~ 12 , n = 9. Figure |2] shows the locations of nodes in [—3,3], 
and the relative error of the quadrature as a function of k G [1/16,4]. 
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Figure 2: Quadrature nodes and relative error for G(k,x) = x k 
4.3 Exponentials k x , hyperbolic Gramians 

This subsection is analogous to the preceding one; therefore, we will only provide the 
essentials. The family of exponential functions 

G(k,x) = { k x , xe[a,b], fee [a,/?]}, a>0 (84) 

is not equivalent to exp(kx). The minimal function is dependent on k but the dependence 
is separable 

fj,(x, k) = x/k := /i(x)/k (85) 

The factor space 

T(k,x) — { k x , xe[a,b], k e [y/a, \f]3\ } (86) 
gives rise to the folding kernel 

F(k,k',K) = 5(kk' -k), k,k' e [Va, y/0\, Ke[a t 0\} (87) 

Therefore, the Gramians A(k, k') and B(k, k') are operators with ks'(k) and s(k) on the 
hyperbolae kk' = k. For constant weight u — 1, 

„, fc fe -A; a , bk b -ak a s(k) 

s ( k > = i I » fcs = — T~T 

In /c In k In fc 

18 



For a numerical experiment, we construct a Gaussian quadrature for G(k,x) = k x , 
x G [a, b] = [—3,3], k G [a, /3] = [1/16,4] by constructing an inner product quadrature 
for the factor space T(k,x) = k x , x G [a, b] = [-3,3], k G [y/a,y/j3\ = [1/4,2]. For a 
cut off precision e = 10~ 12 , the procedures of §3.41 gives rise to n = 9. Figure E] shows 
the locations of nodes in [—3,3], and the relative error of the quadrature as a function 
of k G [1/16,4]. 

Nodes for k x ; xe [-3,3], ke [1/16,4] 




x 10" 



-10 12 
Relative error as function of ke [1/16,4] 




Figure 3: Quadrature nodes and relative error for G(k,x) = k x 



5 Generalizations and applications 

The algorithms for the inner product quadrature design, presented in Theorem 
I2.3[ and I3.5[ will also work for matrix and tensor quadrature weights. Take the two 
Gramians A, B of Theorem I3.5I for example, the product space T1(S, T) of (|35|) may have 
a dimension on the order mn. A quadrature of r nodes, with r < m, will integrate these 
0(mn) distinct functions only if the quadrature weights w has off diagonal entries. It 
may be a banded matrix, or a dense matrix with a predetermined diagonals, but as long 
as the r-by-r matrix w is invertible, Theorem I3.5I still holds, for its proof is equally valid 
as the diagonal matrix w is replace by an invertible one. 
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Tensor "weights" refer to one or two r-by-r matrix w which will entrywise multiply 
the integrand T r from the right, or S r from the left, or both, as opposed to standard 
matrix-matrix multiplication. The proof of Theorem 13.51 holds, as T r and S r will still be 
full rank r after the entrywise multiplication, otherwise the rank of B will be less than 
r. 

Matrix and tensor quadrature weights are related to certain sensing and inverse scat- 
tering applications; see §5.11 for more details. 

The 1-D results presented in this paper makes a step toward a systematic method 
to design Gaussian quadratures for an arbitrary system of functions in one and higher 
dimensions; see §5.21 and 15.31 for 2-D extensions. 



5.1 Separation principle of imaging 

The mathematical models for imaging, with the notable exceptions of MRI and X-ray 
CAT scan due to absence of wave scattering as their probing mechanisms, are inconsistent 
in that their formulation is based on reflectivity or scattering coefficient of targets as a 
function of position. But in many applications, these functions are not nearly single 
valued. Amplitude of backward, monostatic reflected wave from a small target depends 
on direction unless the target is a ball, for example, with uniform reflection coefficient 
on the sphere. 

There is a remarkable property of Gaussian quadrature design - the nodes can be 
determined first and independently of the weights. This is also valid for a "quadrature" 
with inconsistent "quadrature weights", namely with tensor weights. For imaging or 
inverse scattering with waves, the measurement is typically a Gramian matrix known 
as the scattering matrix. For some r point targets as the scatterers, there is a r-term 
quadrature to integrate the Gramian matrix, and the quadrature notes fall on the loca- 
tions of the point targets, provided that the size of the Gramian matrix is no less than 
r. Thus, the quadrature approach presents an alternative model based on the locations 
of targets. 

If we construct a quadrature for the Gramian matrix, the locations of the targets 
will be determined first and separately from the target's reflectivities, whether or not 
they are consistent. If consistent, and if there is no multiple scattering among them then 
the quadrature weights will be the reflectivities; if there is multiple scattering then the 
quadrature weights will be a dense matrix which together with the quadrature nodes 
will be sufficient to recover the consistent reflectivities via solution of a simple matrix 
equation. If the reflectivities are inconsistent, the quadrature weights will be tensor, and 
it is possible to assign an average reflectivity to each point target. 
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5.2 Quadratures in higher dimensions 

A Gaussian quadrature in two dimensions integrating the bivariate polynomials of degree 
less than 2n in a domain D, as is well known, is a summation of n(n + l)/2 terms. Such 
a quadrature rarely exits. We will, however, provide a 2-D versions of Theorems 12.11 
and 12.21 to construct the quadrature by eigen decomposition, and to illustrate what is 
required of quadrature design in higher dimensions. The results will also be useful in §5.31 
for quadrature in two and higher dimensions constructed by a technique called deflation. 

Let T(n, x, y) of size n(n+l) /2-by-D be the n(n+l)/2 basis functions for polynomials 
of degree less than n in the domain D. Let 

B = T(n,x,y) -T(x,y,n) (89) 
A x = T(n,x,y) ■ xT(x,y,n) (90) 
A y = T(n,x,y) -yT(x,y,n) (91) 

where the dot product is over domain D and with a weight function u. We have 

Theorem 5.1 // there is a n(n + l)/2-term quadrature {{xj,yj)]Wj} to integrate the 
Gramian matrices B, A x , and A y , then 

XjiAxB- 1 ) = Xj , j = l:n(n + l)/2 (92) 
XjiAyB' 1 ) = Vj , j = l:n(n+l)/2 (93) 

Here the weight u is not assumed positive. This result will be useful in §5.21 for deflating 
the Gramians. 



Theorem 5.2 Let the weight function u be positive definite, and let Vj denote the j-th 
eigenvector of a matrix. The three conditions are equivalent 

(i) There is a n{n + l)/2-term quadrature {(xj,yj);vjj} to integrate B, A x , and A y . 

(ii) The two quotient matrices share common eigen space, and 

X^B- 1 ) = x v j = l:n(n + l)/2 (94) 
XjiAyB' 1 ) = y j} j = l:n(n + l)/2 (95) 
V] (A X B- 1 ) = Vj(A y B~ l ) — T(n,Xj,yj), j = l:n{n + l)/2 (96) 

(Hi) Then + 1 orthogonal polynomials of degree n have n(n+l)/2 real, pairwise distinct, 
common zeros {(xj,yj), j = l:n(n+ l)/2}. 

The proof of equivalency of (i) and (ii) is similar to that of Theorem 12.21 For (iii), see 
the proof of Theorem 12.51 
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5.3 Deflation for 2-D quadrature design 

A node of a 2-D quadrature provides 3 parameters {{xj,yj);Wj}. Denote by pffi the 
linear space of bivariate polynomials of degree less than n. Therefore, 

dim(Pf } ) = n(n + l)/2, and dim(P 2 ( J) = n(2n + 1) (97) 

A quadrature integrating P 2n generally requires no less than a third as many nodes as 
the dimension, namely n(2n + l)/3 nodes. 

A classical Gaussian quadrature for bivariate polynomials, if exists, can be con- 
structed by Theorem 15 .![ using n(n + l)/2 nodes to integrate P 2n > therefore, the quadra- 
ture problem is over-determined and rarely has a solution. The algorithm of Theorem 
15.11 is rarely useful. But it can be modified and made useful by deflating the Gramians 
A and B iteratively. 

The eigen decomposition of Theorem 15.11 can only provide n(n + l)/2 nodes. Ad- 
ditional nodes will be determined by other mechanisms. The number of these nodes 
is 

dN = n{2n + l)/3 - n(n + l)/2 = n(n - 1) /6 (98) 

which is about a third of n(n + l)/2, namely a third of what can be provided by the 
eigen decomposition. In 3-D, the ratio is 1; as many additional nodes are requires as 
those by the eigen decomposition. Deflation is a method to provide the additional nodes 
iteratively. The following description takes bivariate polynomials in a triangle as example 
to illustrate the method. 

1. Suppose that a total of 40 nodes are required to integrate polynomials of degree less 
than 2n for some n. Suppose that the size of Gramians is 30, so eigen decomposition 
can only provide 30 nodes. Additional 10 nodes will be supplied by an iterative 
procedure. 

2. Suppose we are given the precise locations of 10 out of the 40 nodes and the 
corresponding weights Wj. Each node Zj = (xj, i/j), j = 1 : 10, gives rise to a rank 
one matrix T(n, Zj)wjT(zj,n); see Theorem 15.11 for notation. Deflation involves 
three steps (i) Remove these 10 matrices from Gramian B (ii) Remove the 10 rank 
one matrices T(n, Zj)wjXjT(zj,n) from Gramian A x (iii) Remove the 10 rank one 
matrices T(n, Zj)wji/jT(zj,n) from Gramian A y . 

3. Theorem 15.31 below states that the eigen decomposition on the quotient matrices 
A X B~ X and A y B~ x (cf Theorem 15. ip after the deflations will provide the exact 
locations of the remaining 30 nodes. 

4. Initialization. Choose 10 nodes and weights as initial guess. There are ways to make 
good initial guess located in a corner of the triangle - the domain of integration. 
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5. Iteration. Eigen decomposition of the deflated quotient matrices to obtain 30 nodes. 
Discard 20 of them by choosing only 10 out of 30 that are farthest from the 10 
initial guess, and use them as the initial guess for the next iteration. 

6. Convergence. The Coulomb potential 1/r decays over distance. Its perturbation 
due to that of charge location decays faster: 1/r 2 . The location errors in the 10 
initial guess will have minimal influence on the farthest of the 30 nodes. 

Deflation is also useful for constructing (i) Gauss- Radau type formula (with an end x=a 
or b fixed as a quadrature node) in one and higher dimensions (ii) Gauss-Lobatto type 
formula (with two ends fixed as quadrature nodes) in one and higher dimensions. 

Deflation can be used for constructing a Gaussian quadrature in a submain and 
merging it to an existing quadrature as the trapezoidal rule in another subdomainsuch - 
the so-called hybrid rules [3]. 

Theorem 5.3 Suppose there is a n(n + l)/2 + r-term quadrature {(xj,yj);Wj} to inte- 
grate the Gramians B, A x , and A y of / fgPj) -/ fP7]) . For the first r nodes, let the deflated 
Gramians be defined by 

r 

B = B - ^T(n,x i ,y i )w; i T(x i ,y i ,ra) 

r 

A x = A x - T(n, Xj, y j )w j x j T(x j ,y j , n) 

3=1 
r 

3=1 

then 

XjiAxB- 1 ) = Xj , j — 1 + r:n(n + l)/2 + r (102) 

XjiAyB- 1 ) = Vj , j = 1 + r:n(n + l)/2 + r (103) 

The proof is a direct consequence of Theorem 15.11 applied to the deflated weight function 

r 

u(x) = u(x) - ^ w 3 5 ( x - x jiV- Vj) ( 104 ) 

3=1 
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